Critical temperature and the transition from quantum to classical order parameter 
fluctuations in the three-dimensional Heisenberg antiferromagnet 
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We present results of extensive quantum Monte Carlo simulations of the three-dimensional (3D) 
5=1/2 Heisenberg antiferromagnet. Finite-size scaling of the spin stiffness and the sublattice 
magnetization gives the critical temperature T c / J = 0.946±0.001. The critical behavior is consistent 
with the classical 3D Heisenberg universality class, as expected. We discuss the general nature of 
the transition from quantum mechanical to classical (thermal) order parameter fluctuations at a 
continuous T c > phase transition. 
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The Heisenberg quantum antiferromagnet [lj,U is one 
of the most important models in solid state physics. A 
wide variety of physical behavior can be realized, depend- 
ing on the geometry of the lattice over which the pair 
interaction JSi ■ Sj is summed. On a ID linear chain, 
the quantum fluctuations destroy the long range order, 
leading to a ground state with gapped (for integer S) or 
critical (for half odd-integer S) fluctuations ||-||] • On a 
2D square lattice the ground state is ordered for all S 
HQ . For T > the order is destroyed by thermal fluctu- 
ations ||. In 3D, a second-order phase transition to an 
ordered state at T c > is expected. The quantum fluctu- 
ations of the order parameter become negligible (relative 
to the classical, thermal fluctuations) in the neighbor- 
hood of T c and the critical behavior should therefore be 
that of the classical 3D Heisenberg model. The quan- 
tum fluctuations vanish also in the limit T — > oo, and 
hence their relative strength should have a maximum for 
T c < T < oo. A universal "quantum critical" regime 
|p[-|lT| is in principle possible if the system is sufficiently 
close to a quantum critical (T c = 0) point jL2j. 

Analytic calculations, employing several different ap- 
proximations at finite T, have recently been carried out 
for the S = 1/2 Heisenberg model on a 3D simple 
cubic lattice Critical temperatures in the range 

T c / J = 0.89 — 1.13 were obtained. Results for thermody- 
namic quantities depend strongly on the approximation 
used. High temperature series expansions have been used 
for the thermodynamics at high temperatures |1J], but 
cannot be expected to be accurate close to T c . Quan- 
tum Monte Carlo (QMC) methods have been important 
for non-approximate calculations for 2D quantum anti- 
ferromagnets, but so far only limited results have been 
obtained for 3D systems fTq |. We here present detailed 
QMC calculations aimed at reliably extracting T c as well 
as providing unbiased results for the temperature depen- 
dence of thermodynamic quantities. We also discuss the 
precise nature of the crossover from quantum mechanical 
to purely classical thermal order parameter fluctuations 
at the critical point. In particular, we show that the 
ratio 5(Q)/[Tx(Q)] of the static order parameter struc- 



ture factor and susceptibility obeys a finite-size scaling 
at T c that is completely universal, independent of the ex- 
ponents of the phase transition or the dimensionality. 

For the QMC simulations we have used the Stochastic 
Series Expansion (SSE) algorithm [jl6| (a generalization 
of Handscomb's method p7j ) for L 3 lattices with even L 
up to 16 in the grand canonical ensemble. The method 
does not introduce any systematic errors in computed 
quantities. As a significant development allowing us to 
efficiently obtain results on a dense temperature grid, we 
have introduced "tempering" within the SSE algorithm. 
In this approach [jl8j , the configuration space of the sim- 
ulation is extended to include a range of temperatures 
between which the system can make transitions, hence 
allowing calculations for a large number of temperatures 
much faster than with separate fixed-T simulations. Also 
crucial in this work is the use of covariance effects to in- 
crease the numerical precision. The scheme introduced 
in |l{J allows for reduction of the statistical errors by al- 
most three orders of magnitude for the 3D Heisenberg 
model. 

In the path integral formulation, a quantum system 
at finite T corresponds to a classical system with an ex- 
tra "imaginary time" dimension of length L T ~ 1 /T [Q . 
Standard finite-size scaling arguments plj can be applied 
to this system. In order to determine the critical temper- 
ature T c , we here first consider the T and L dependence 
of the spin stiffness p s . With the SSE method p s can 
be calculated directly from the global winding number 
fluctuations in periodic systems , as described in Q . 
Hyperscaling theory predicts the finite-size scaling of the 
stiffness at the critical point yM, 



p s {T c ) = L 



2 — rf- 



(1) 



where d is the dimensionality and z is the dynamic critical 
exponent which is zero for a T c > transition. For the 
3D Heisenberg model, Lp s graphed versus T for different 
L should hence intersect at T c . For L = 4 — 16 the 
curves indeed intersect almost at a single point, as shown 
in Figure For the smallest systems there is a clear 
tendency of the intersection points for sizes L and L + 2 



1 



to move to lower T with increasing L, but for L > 10 
no shifts can be seen within the statistical errors. From 
these results we therefore estimate T c w 0.945 J. 
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FIG. 1. Spin stiffness vs. temperature for systems with 
L = 4 — 16 (the slope increases with L). 

In the thermodynamic limit, the sublattice magnetiza- 
tion |(<Sf)| (the order parameter) increases continuously 
from zero for T < T c . Consider the static structure factor 
and susceptibility 



%)4E e<q ' (rrri) ( s i s ^ 



(2a) 



1,3 



^) = ^E eiq ' (rj " rs) / dT(S?(r)S?(0))- (2b) 

At the ordering wave- vector Q = (tt, it, 7r), S and \ 
should both scale with the system size as L 2 ~ v at T = T c , 
where r\ w 0.03 for the 3D Heisenberg transition p5| . 
%(Q) corresponds to the full space-time integral of the 
sublattice magnetization of the 3+1 dimensional classical 
system and can therefore be expected to be less affected 
by corrections to scaling. Figure ||a shows In (x/L 2 ) 
vs. In (L) at temperatures close to the T c estimated above 
using the stiffness. Asymptotically we expect the points 
to fall on a straight line with slope —77 w —0.03 if T = T c , 
and eventually diverge upward or downward below and 
above T c , respectively. At T w 0.946 J the L > 10 data 
show the expected scaling. Investigating the behavior 
for several temperatures in this neighborhood, and us- 
ing the estimated accuracy of the exponent 77 fl25fl , gives 
T c /J = 0.946 ± 0.001, in excellent agreement with the 
spin stiffness scaling. 

With T c determined, the expected scaling for T > T c 
can be tested. In the thermodynamic limit x should di- 
verge as i -7 , where t = \T — T c \ and 7 = v{2 — rj) w 1.40 
Finite-size scaling predicts Xh{t) = Xoa(t)f[£(t)/L], 
where the correlation length £ ~ t~ v . Hence, xz,(t)t 7 
graphed versus Lt" for different L should collapse onto a 
single curve. As shown in Figure ||b, the data are indeed 
collapsed, confirming the expected universality class and 
the estimated T c . 
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FIG. 2. (a) Size dependence of the staggered susceptibil- 
ity divided by L 2 at temperatures close to T c . The straight 
line is the L — > 00 behavior (slope —0.03) expected with the 
exponent 77 = 0.03. (b) Finite-size scaling of the staggered 
susceptibility above T c , using the classical 3D Heisenberg ex- 
ponents 77 = 0.03 and v = 0.71, and assuming T c /J = 0.9459. 

We next consider the specific heat, C — (dE /&T)/ L 3 . 
Figure [s|a shows the temperature dependence for different 
system sizes. The position of the maximum of C repre- 
sents a size dependent definition of the critical point, and, 
as shown in more detail in Fig. ||b, seems to indicate that 
T c is lower than the value estimated above. We believe 
that the reason for this apparent inconsistency is to be 
found in the fact that changing T of the quantum system 
corresponds to changing L T of the corresponding 3+1 di- 
mensional classical system. The calculated C hence cor- 
responds to dE/dL^ 1 . The critical behavior, governed 
by the exponent a, of the actual classical specific heat 
and C of the quantum system should nevertheless be the 
same (except exactly at a T c = quantum critical point) 
since the boundary between the ordered and disordered 
phases in the classical (T, L T ) plane is analytic for T > 0. 
However, for small systems the shift of the peak position 
with L is relatively large, and the associated L T variation 
can then cause significant deviations from the asymptotic 
finite-size scaling behavior of C. In order to partially 
compensate for the volume change due to the variable 
L Tl we consider the corresponding quantity defined per 
full space-time volume, i.e., TC. Figure yc shows the 
behavior of the peak of TC. Compared to the behav- 
ior of C, we note that the peaks for these small systems 
are shifted noticeably to the right, and now appear to be 
consistent with T c k, 0.946 J. 
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FIG. 3. (a) Specific heat vs. temperature for L = 4 — 16. 
The peak height increases with L. (b) the behavior of the 
peak for L — 6 — 16 on a more detailed scale, (c) The peak 
behavior of TC. The estimated T c is between the dashed 
lines. 

For the 3D classical Heisenberg transition a is negative 
(a w —0.11) [ p5[ ], and hence C should exhibit a cusp- 
like singularity maximum as L — > oo. Finite-size scaling 
predicts the size dependence of the maximum value as 
C max (L) = C max (oo) - aL a l v (v « 0.71). The data for 
TC is indeed well described by this expression for L > 6, 
and gives T c C max (oo) ea 2.4. The finite-size shift of the 
peak temperature should scale as L~ x l v . The results do 
not show this behavior. It is quite likely that this is 
due to the L T effect discussed above, which is not fully 
eliminated in TC. We note, however, that the expected 
scaling of C is not easily observed even for the classical 
3D Heisenberg model, for which significantly larger sys- 
tems have been studied [25j| . A qualitative difference is 
that the peak shifts to higher T with increasing L in the 
classical case [p5[ . As discussed above, for the quantum 
model the peak in C initially moves with increasing L 
from above T c to slightly below T c . As L — > oo, C and 
TC should both peak exactly at T c . One can therefore 
expect a change in the sign of the C peak shift for some 



large L, and then the same asymptotic behavior as for 
the classical model. 

We now discuss the nature of the crossover to purely 
classical (thermal) fluctuations of the order parameter as 
T — > T c from above. For this purpose we note the well 
known sumrulcs 



1 f 00 

S(q) = - / dw(l + e 

7T Jo 

2 f°° 

X (q) = - / duu-^l 
t Jo 



)S(q,w), 



-f3uj 



)S(q,w), 



(3) 
(4) 



where 5(q, to) is the dynamic structure factor, and define 
the ratio 



i?(q)=5(q)/[T X (q)]. 



(5) 



For purely thermal fluctuations of momentum q, all the 
spectral weight is in a 6 function at u> = 0, and the clas- 
sical relation -R(q) = 1 is obeyed. Quantum (\u>\ > 0) 
fluctuations increase R above 1. For the 3D antiferro- 
magnet the fluctuations at the antiferromagnetic momen- 
tum Q should be purely classical for T < T c . As T — > oo, 
i?(q) — ► 1 for any q. For T T c from above, L T <C £ and 
in evaluating Eqs. (^) the imaginary-time dependence of 
the long-distance correlation function can then be ex- 
panded to leading order in r/r (or r 1 / 2 fr if the system 
is not Lorenz invariant, which does not change the re- 



sult). This gives R(Q) - 1 ~ £~ 



t , and the size 



dependence at T 



T c should therefore be (substituting 



R(Q) — 1 ~ L" 



(6) 



Note that this scaling does not contain any model- 
dependent exponents, and remains valid in any number 
of dimensions (with the S'(Q) and x(Q) corresponding 
to the relevant order parameter) as long as T c > and 
the order parameter does not commute with the Hamil- 
tonian. (In the commuting case i?(Q) = 1 for all L and 
T). R(Q) should therefore in principle be useful, as an 
alternative to the stiffness scaling (|l|), for extracting T c 
in cases where the critical exponents are not known. 

Figure |a shows R(Q) - 1 versus temperature for 3D 
Heisenberg lattices. The finite-size gap leads to a T — > 
divergence starting at a temperature which decreases 
with increasing L. For the larger systems there is a max- 
imum at T w 1.7 J, and a rapid drop below this tem- 



perature. The maximum value R n 



1 « 0.05, which 



is about half of the value for the 2D Heisenberg model 
Jl0| , pr|| , reflecting the reduction of quantum fluctuations 
with increasing dimensionality. 

Figure. ||b shows the T dependence of [-R(Q) — l]L 2 
for different system sizes. The curves indeed intersect. 
The intersection between L and L + 2 moves to lower T 
with increasing L, and is consistent with the estimated 
T c w 0.946 J. The proposed scaling law, Eq. @, for R(Q) 
is hence confirmed. 
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FIG. 4. (a) Temperature dependence of R(Q) — 1 for system 
sizes L = 4 — 16 (values decrease with increasing L). (b) The 
same quantity multiplied by L? close to T c (slope increases 
with L). 

It is well known that a state of dominant classical order 
parameter fluctuations can extend to a finite regime well 
above T c §JT|. This is not in conflict with the universal 
scaling behavior of i?(Q) noted here, according to which 
the classical value i?(Q) = 1 always is approached alge- 
braically as T — > T c . To see this, consider a system of 2D 
layers coupled by a very weak interlayer coupling J± . As 
the temperature is lowered, the correlation length within 
the layers initially grows exponentially || but there are 
no significant correlations between the layers. For T suf- 
ficiently low, but T>T C , they are in the "renormalizcd 
classical" Q regime and R(Q) — > 1 exponentially [ p6[ . 
However, as T is further lowered to the point, close to 
T c , where the layers start to correlate with each other, 
there will be a crossover to a £ ~ t~ v scaling behavior 
and the form (§) for R(Q). At this point R(Q) is al- 
ready exponentially close to 1, and the further algebraic 
reduction of R(Q) as T — ► T c is exponentially small. 

In summary, we have determined the critical tempera- 
ture T c /J = 0.946±0.001 of the 3D Heisenberg antiferro- 
magnet with S = 1/2, using QMC data of high precision. 
We have discussed the general nature of the crossover to 
purely classical thermal order parameter fluctuations as 
T — > T c from above, and shown that asymptotically the 
transition is universal in the ratio S(Q)/[Tx(Q)], inde- 
pendent of the universality class of the critical point. We 
have also noted how quantum mechanical effects can af- 
fect the finite-size scaling behavior of the specific heat. 
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